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Abstract In order to design a WSB trajectory to Moon, 
highly elliptical geocentric orbit is propagated forward in 
time so that its perigee increases to Earth-Moon distance. 
Another highly elliptical selenocentric orbit is propagated 
backward in time till it starts moving towards Earth. The 
two trajectories are patched using Fixed Time of Arrival tar¬ 
geting method to obtain a WSB transfer trajectory. In this 
paper, various geocentric and selenocentric orbits are con¬ 
sidered and eligible candidates for WSB transfer are repre¬ 
sented on the phase space with color code on transfer time. 
Genetic Algorithm is used to find patching points to reduce 
the sum of AV required at the patching points. Lunar fly-by 
on the way to apogee is also studied. 

Keywords Weak stability boundary • Fly-by • Fixed time 
of arrival targeting • Genetic Algorithm 

1 Introduction 

Belbruno (1987) discovered a new method to design low en¬ 
ergy Earth-to-Moon transfers using weak stability boundary 
(WSB) at Moon. WSB is region in the phase space where the 
perturbative effects of Earth-Moon-Sun acting on the space¬ 
craft tend to balance. These transfers take the advantage of 
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Sun’s gravity to increase perigee altitude of a highly ellip¬ 
tical geocentric orbit, from 200-2000 km to Earth-Moon 
distance (384,400 km). The spacecraft reaches Moon with 
sufficient energy to get temporarily ballistically captured by 
Moon. Specific Earth-Moon-Sun geometric configuration is 
important to accomplish such transfers. The AV savings is 
of the order of 18 % compared to conventional Hohmann 
transfer and time of flight is between 60 to 100 days (Bel¬ 
bruno and Miller 1993). Japanese Hiten Mission in 1991 
successfully demonstrated WSB transfer to Moon. 

The study of restricted three- and four-body problems 
(Belbruno 2004; Koon et al. 2006; Parker and Anderson 
2014) have been used, in recent years, to design low en¬ 
ergy transfer trajectories or WSB transfer trajectories. These 
highly non-linear orbits obtained due to dynamics of 3 
bodies require less propellant compared to conventional 
patched-conic orbits but transfer time is longer compared 
to latter. Other advantages include extended launch period 
and gradual capture or relaxed operational schedule. On the 
other hand, these transfers require high reliability of sub¬ 
systems and long link distance. For these reasons these tra¬ 
jectories are mostly suitable for rescue mission, cargo mis¬ 
sions for permanent lunar base, sample return etc. and not 
suitable for manned missions. Hiten in 1991 (Belbruno and 
Miller 1990; Uesugi 1991), SMART-1 in 2003 (Foing and 
Racca 1999), ARTEMIS in 2010 (Broschart et al. 2009) and 
GRAIL in 2011 (Roncoli and Fujii 2010; Chung et al. 2010; 
Hatch et al. 2010) have traversed low-energy transfers from 
Earth to Moon. 

Since their discovery in 1987, WSB transfers were ex¬ 
plored further by Belbruno (1990, 2004, 2008), Belbruno 
and Miller (1993), Miller and Belbruno (1991), Yamakawa 
(1992), Yamakawa et al. (1993), Koon et al. (2000, 2006), 
Ivashkin (2002, 2003, 2004), Parker and Lo (2005), Parker 
(2007), Parker and Born (2008), Parker (2010), Parker et al. 
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(2011), Parker and Anderson (2014), Topputo et al. (2005), 
Topputo (2013) and many more. 

Yamakawa (1992) describes Gravitational Capture as 
“Although local C3 with respect to Moon is positive out¬ 
side the sphere of influence, finally it becomes nega¬ 
tive at perilune”. The perilunes located in first and third 
quadrant in selenocentric Earth-Moon fixed frame can be 
achieved by gravitational capture from a geocentric orbit. 
Solar gravity helps to enlarge the perigee from LEO dis¬ 
tance (200-2000 km) to Earth-Moon distance and to at¬ 
tain Moon from large geocentric orbit with semi-major axis 
over 500,000 km. Spacecraft located in second and fourth 
quadrant in Sun-Earth fixed frame yields increase in local 
perigee distance. Lunar swing-by reduces local eccentricity 
(and helps in perigee raise) by which total flight time for 
Earth-Moon transfer can be reduced. It also gives flexibil¬ 
ity in orbit design by controlling initial semi-major axis and 
swing-by distance. 

In this paper, lunar capture trajectories are obtained by 
back-propagation in the restricted three-body (Earth-Moon- 
spacecraft) problem and are represented in phase-space with 
color code on time of capture. Highly elliptical geocentric 
orbits are investigated in Sun-Earth-spacecraft system and 
those for which the perigee increases to Earth-Moon dis¬ 
tance are also represented in phase space with “color code 
on time of flight”. Time of flight is the orbital period of 
the elliptical geocentric orbits, the time required to raise the 
perigee altitude to Earth-Moon distance. By “color code on 
time of capture” we mean that different colors have been as¬ 
signed to different ranges of time of capture (number of days 
the trajectory from an orbit around moon is back propagated 
so that its C3 wrt moon becomes positive), so that looking 
at the figure short flight duration trajectories can be differen¬ 
tiated from longer flight duration trajectories. Moreover the 
distribution of trajectories with different capture durations 
in the phase space can be clearly visualized. The two trajec¬ 
tory segments are joined using Fixed Time of Arrival Target¬ 
ing (FTAT) method. Genetic Algorithm (GA) is used to find 
patching points to reduce AV. The role of initial phase angle 
for perigee enlargement is highlighted. Also lunar swing-by 
on the way to apogee is studied. 

This study will be helpful for mission planners to design 
a WSB trajectory of their choice of arrival point and flight 
duration. Due to color code on time of capture in the phase 
space diagrams, an approximation of total flight duration can 
be made without actually constructing the complete trajec¬ 
tory. This method is different from other methods available 
in literature because it provides feasible options for arrival 
(lunar capture orbit) and departure (highly elliptical geocen¬ 
tric orbit) conditions for the construction of a WSB trajec¬ 
tory from Earth to Moon. Given departure and arrival condi¬ 
tions, GA is used to find patching points on these trajectories 
so that the patching AV is minimized. Moreover due to the 
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detailed study of arrival and departure conditions, the algo¬ 
rithm to design WSB trajectory, presented here, is designed 
in such a way that the dimension of search space is reduced 
to two. 


2 Equations of motion 


Lunar capture trajectories are obtained in Circular Restricted 
Three-Body Problem (CR3BP) by back propagation of ini¬ 
tial conditions near Moon. The equations of motion of 
CR3BP are given in Appendix. Fixed time of arrival target¬ 
ing, as given in Appendix, is implemented to patch two tra¬ 
jectories, namely, one forward propagated from Earth and 
other capture trajectory obtained by back-propagation. Ge¬ 
netic algorithm is implemented to find patching points on 
the two trajectories. During forward propagation from Earth, 
the role of Sun to raise perigee altitude is investigated in Re¬ 
stricted Four Body Problem (R4BP). The equations of mo¬ 
tion of R4BP are given in Appendix. Also fly-by Moon on 
the way to WSB Earth to reduce time of flight is evaluated. 

The position and velocity of a particle in orbit around 
smaller primary m 2 are given below (Garcia and Gomez 
2007). 

Initial conditions with positive velocity (osculating retro¬ 
grade motions about m 2 ) 

x = — 1 + /z + r 2 cos 0 , y = 7*2 sin 0, (1) 

x = r 2 sin0 — v sin0, y = — 7*2 cos 0 + vcosO. (2) 

Initial conditions with negative velocity (osculating direct 
motions about m 2 ) 

x = — 1 + fi + 7*2 cos 0, y = 7*2 sin 0, (3) 

x = 7*2 sin0 + usin0, y = — 7*2 cos 0 — vcosO, (4) 


where 7*2 is the distance between m 2 and m3, 0 is the angle 
between radial line and x-axis and /z is the mass ratio. In 
both cases, the modulus of the initial sidereal velocity of the 
infinitesimal particle is equal to v. The value of v is given 
by 



/z(l + e) 


( 5 ) 


The above equations are used in Sect. 3.1 to obtain initial 
conditions for back propagation to find lunar capture trajec¬ 
tories. 

In order to understand the role of Sun to enlarge perigee 
altitude, initial conditions near Earth are propagated forward 
in time in the R4BP. The two body orbital elements in planar 
case associated to the orbit of m3 around mi, are given by 
the following map (Belbruno 2008) 


(a,e,co,(p) -> (x,y,x,y). 


( 6 ) 
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Fig. 1 Lunar capture 
trajectories obtained by back 
propagation for 50 days. The 
initial conditions which lead to a 
capture trajectory are plotted in 
the phase space with color code 
on time of capture, (a) Direct 
motion about m 2 , (b) retrograde 
motion about m 2 
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Here a is the semi-major axis, e is the eccentricity, co (re¬ 
ferred as AOP later in the text/figures) is the angle between 
v-axis and line of apsis and ep is the angle between line 
of apsis and radial line in m \ centered co-ordinate system. 
At t = to, when the rotating and the inertial Earth-centered 


frames are parallel, we can write 

x = r cos(&> + ep) — /x, (7) 

y = r sin(&> ± ep), (8) 

x =r cos(&> + ep) — rep sin(&> + <p) + r sin(&> + ep), (9) 

y = r sin(&> ± (p) + rep cos (co ± ep) — r cos(&> ± ep), (10) 

where 


i.e., the time of back-propagation till C3 wrt Moon becomes 
positive, is recorded in each case. Figure 1(a) and (b) shows 
the phase space representation of lunar capture trajectories 
with color code on time of capture. Similar trend is observed 
with perilune 500 km. With the help of these phase space 
diagrams, one can select an appropriate arrival point with 
respect to desired time of flight. The trend observed for di¬ 
rect motion is well within —55° to 55° and 125° to 235° 
AOP as predicted by Yamakawa (1992), Preposition 3.3. In 
case of retrograde motion for higher apolune altitudes the 
above mentioned trend is violated to some extent. This may 
be due to truncation of Taylor’s series during the derivation 
of Preposition 3.3. 


a(l-e 2 ) 

r =-, 

1 -he cos ep 

y/a{ 1 - e 2 )(l - fi) 
<P = - 2 - 

ae(\ — e 2 )epsmep 

r =-. 

(1 + ecos ep) 2 


( 11 ) 

( 12 ) 

(13) 


3 Dynamics of WSB trajectories 
3.1 Lunar capture trajectories 

The model used here is Restricted Three Body Problem 
(R3BP)—Earth-Moon-spacecraft, to back propagate initial 
conditions obtained from (l)-(5) and using equations of mo¬ 
tion provided in Appendix. Perilune is fixed at 100 km and 
apolune is varied from 100 km to 50,000 km. The argument 
of perilune is varied from 0° to 360°. Perilune of 100 km is 
selected as it gives the most useful mapping orbit around 
Moon. Most of lunar missions including Chandrayaan-1 
have been placed in 100 km circular orbit around Moon. 
Apolune is varied till 50,000 so that the eccentricity of the 
100 x 50,000 km orbit is 0.93 which is approximately where 
ballistic capture occurs. Once captured the orbit apolune can 
slowly be decreased to get a circular orbit. Time of capture, 


3.2 Enlargement of perigee from LEO to Earth-Moon 
distance 

Using (6)-(13), Earth Parking Orbits (EPO) with perigee al¬ 
titude 200 km, apogee altitude varying from 1 to 1.5 mil¬ 
lion km (distance of LI point from Earth in Sun-Earth Sys¬ 
tem) and argument of perigee varying from 0° to 360° are 
obtained and propagated in Restricted four Body Problem 
(R4BP)—Sun-Earth-Moon-S/C, for one revolution using 
equations given in Appendix. The perigee altitude is fixed 
at 200 km as most launch vehicles inject their payload in El¬ 
liptical Earth Parking Orbits with perigee at about 200 km. 

The orbits for which the perigee altitude increased from 
200 km to 384,400 ± 50,000 km were represented in the 
phase space in Fig. 2, with color code on time of flight, along 
with some orbits. ±50,000 km is taken so that the trajectory 
is within the sphere of influence of Moon. The initial lu¬ 
nar phase angle ($mo) is 0°. In the phase space diagram two 
peaks are visible in the 2nd and 4th quadrant. This supports 
the result by Yamakawa (1992) that S/C location in 2nd and 
4th quadrant in Sun-Earth fixed frame yields increase in lo¬ 
cal perigee distance. This diagram is also useful for mission 
designers as one can identify appropriate departure condi¬ 
tion with respect to time of flight. 

There is a slight variation in the phase space diagram 
when the initial phase of Moon changes. Figure 3 compares 
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Fig. 2 EPO whose perigee rose 
from 200 km to about 
384,400 km propagated in R4BP 
(Initial lunar phase angle is 0°) 
with color code on time of flight 



xIO 6 



Fig. 3 Comparison of phase space representation of EPO whose 
perigee rose from 200 km to 384,400 km for initial lunar phase angle 
(<9 M0 ) 0° and 90° 

the phase space representation of EPO whose perigee rose 
from 200 km to about 384,400 km for initial phase angle 0° 
and 90°. This tells us that the choice of initial lunar phase 
angle is important for WSB trajectory design. Figure 4 indi¬ 
cates the EPO which encounter Moon when the s/c returns 
to the perigee. All the EPO for which the distance between 
Moon and s/c at perigee is less than 50,000 km are plotted. 
Here the initial lunar phase angle considered is 0°. 

3.3 Fly-by Moon on the way to apogee 

The initial lunar phase angle is adjusted so that fly-by oc¬ 
curs on the way to apogee. Figure 5 shows the phase space 


representation of the EPO, with perigee fixed at 200 km, 
apogee varying from 0.4 to 3 million km and AOP vary¬ 
ing from 0° to 360° in each case, which have a fly-by on 
the way to apogee and their perigee increases from 200 km 
to 384,400 ± 50,000 km. It is observed that the opportu¬ 
nities have increased and flight duration has decreased in 
many cases compared to Figs. 2-4. Figure 6 gives the de¬ 
pendence of time of flight on the angular position of Moon 
during flyby on the way to apogee. 


4 Construction of WSB trajectory 

The phase space diagrams provided in Section 3 helps in lo¬ 
cating feasible departure (highly eccentric geocentric orbit) 
and arrival (lunar capture orbit) conditions. Once the depar¬ 
ture and arrival orbits are identified they are patched using 
Fixed Time of Arrival Targeting (FTAT) algorithm (Inner 
Foop) provided in Appendix. 

As an example, Fig. 7 gives Earth-centered orbit that 
fly-by Moon. The initial conditions are: perigee altitude is 
200 km, eccentricity is 0.976, AOP and true anomaly are 0°. 
These initial conditions are propagated forward in time till 
a lunar fly-by is obtained. It takes about 68.7 days to reach 
near Moon. A point A is selected on this trajectory, which 
will be the patching point for FTAT. Another set of initial 
conditions near Moon are propagated backward in time till 
it reaches about 10 5 km away from Moon or its energy 
wrt Moon becomes positive. One such trajectory is shown 
in Fig. 8. The initial conditions (ICs) for this trajectory is 
perilune is 100 km, eccentricity is 0.917 and AOP = 0°. 
It requires approximately 39 days to escape from Moon. 
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Fig. 4 EPO which encounter 
Moon on return to perigee 
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Fig. 5 Fly-by Moon on the way 
to apogee with color code on 
time of flight 



A point B is selected on this trajectory for patching with the 
previous one. At present these points A and B are selected 
judicially and are not optimized for minimum patching AV. 
The points A and B are patched using FTAT. The resultant 
trajectory is a ballistic capture trajectory at Moon as shown 
in Fig. 9, and so theoretically no impulse is required for orbit 
insertion. 

Real-valued Genetic Algorithm (GA) is implemented to 
find patching points to minimize patching AV requirements 
(Deb 1995, 2001). GA provides sufficiently good solutions 
in spite of highly nonlinear dynamics involved in the sys¬ 
tem. 


Variables: 

1. Time t\ is the propagation time from IC near Earth to the 
patching point A. 

2. Time t 2 is the back propagation time from IC near Moon 
to the patching point B. 

Objective is to minimize AV (sum of AV required at both 
patching points A and B). The fitness function is 

1 

(l+Avy 

Characteristics/Parameters of GA: 

1. No. of individuals in a population = 50. 
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Fig. 6a: Fly-by Moon (0Mo-ldeg) 



Fig. 6c: Fly-by Moon (0Mo-O.5deg) 



Fig. 6e: Fly-by Moon (0mo) 



Fig. 6b: Fly-by Moon (0Mo-ldeg) and arrive Moon 
at first perigee passage 



Fig. 6d: Fly-by Moon (0Mo-O.5deg) and arrive 
Moon at first perigee passage 



Fig. 6f: Fly-by Moon (0mo) and arrive Moon at first 
perigee passage 


Fig. 6 Dependence of time of flight on the angular position of Moon during flyby on the way to apogee 


2. No. of generations over which solution evolves = 50. 

3. Crossover probability = 0.85. 

4. Elitism, retaining the best individual in a generation un¬ 
changed in the next generations, is used. 

5. Mutation rate is dynamically adjusted. 

6. Roulette wheel selection is used for crossover. 

In the previous example GA was used to find patching 
points A and B to reduce AV. Figure 10 gives the optimized 
trajectory. 


The algorithm for construction of WSB trajectories using 

the dynamics is as follows: 

1. Select a suitable arrival condition based on the orbit and 
capture duration from Fig. 1. 

2. Select a suitable departure orbit based on orbit, transfer 
time and requirement of lunar fly-by from Figs. 2-6. 

3. The above two segments are joined using FTAT. The 
patching points on the two segments are searched using 
GA. 
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Fig. 6g: Fly-by Moon (0Mo+O.5deg) 



Fig. 6i: Fly-by Moon (0Mo+ldeg) 
Fig. 6 ( Continued ) 



x 


Fig. 7 Earth centered orbit 

Any optimization algorithm can be used instead of GA. 
As a case study, Pattern Search (PS) (Abramson 2002; Au- 
det and Dennis 2003) and Nelder-Mead simplex algorithm 
(Lagarias et al. 1998) were implemented and tested for 6 
test cases. In all these cases, the initial orbit around Earth is 
fixed with perigee altitude 200 km, eccentricity 0.976 and 
AOP and true anomaly are 0°. The perilune altitude is fixed 
as 100 km. For the first three cases initial lunar orbit has 



Fig. 6h: Fly-by Moon (0Mo+0.5deg) and arrive 
Moon at first perigee passage 



Fig. 6j: Fly-by Moon (0Mo+ldeg) and arrive Moon 
at first perigee passage 



Fig. 8 Lunar capture orbit 

AOP 0° and eccentricities are 0.917, 0.928 and 0.938, re¬ 
spectively. For next three cases, AOP of initial lunar orbit is 
180° and eccentricities are 0.921, 0.931 and 0.941, respec¬ 
tively. Table 1 gives the results obtained for these 6 cases 
using three optimization techniques. It is observed that GA 
gives lower total A V for all the cases compared to other two 
methods, while number of function evaluations are lesser for 
PS compared to other two algorithms. 
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Table 1 Comparison between 
optimization methods—pattern 
search, Nelder-Mead simplex 
algorithm and genetic algorithm 


Method 

Case No. 

h (days) 

t 2 (days) 

Total 

AV (km/s) 

No. 

function 

evaluations 

Pattern 

1 

8.97 

16.95 

0.906 

132 

search 

2 

18.33 

20.83 

1.109 

217 


3 

19.57 

42.39 

1.233 

115 


4 

6.81 

28.89 

0.572 

224 


5 

12.34 

0.43 

0.874 

210 


6 

13.13 

16.85 

0.689 

120 

Nelder-Mead 

1 

3.86 

16.93 

0.925 

402 

simplex method 

2 

4.57 

17.40 

1.066 

401 


3 

3.79 

14.27 

1.214 

233 


4 

4.26 

18.65 

1.315 

195 


5 

4.29 

18.92 

1.212 

219 


6 

4.44 

17.57 

1.032 

400 

Genetic 

1 

0.20 

8.33 

0.825 

2600 

Algorithm 

2 

26.15 

8.13 

0.120 

2600 


3 

3.30 

44.99 

0.394 

2600 


4 

43.97 

16.52 

0.418 

2600 


5 

34.90 

5.13 

0.132 

2600 


6 

35.16 

19.46 

0.075 

2600 



Fig. 9 Patching Earth centered orbit with lunar capture orbit using 
FTAT 

Topputo et al. (2004) have presented a method using in¬ 
variant manifold theory to compute transfers from Earth to 
Moon. They start from an Earth Parking Orbit and compute 
trajectory arcs that target a point on using Lambert’s 
three-body problem. A first maneuver places the spacecraft 
into a translunar trajectory starting from an Earth Parking 
Orbit (LEO/GTO); second maneuver injects the spacecraft 
on the capture trajectory on W^i. The total cost of transfer 
is sum of two maneuvers. In an example, they consider two 
departure orbits, 200 km circular LEO and 200 x 35,840 km 
GTO. The given trajectory arrives Moon in an unstable orbit 
around Moon with mean altitude of 21,600 km. The same 



Fig. 10 Patching points A and B obtained by GA 


example has been worked out with the present method us¬ 
ing GA optimization (Table 2). In the present analysis we 
consider an EPO with perigee altitude 200 km, eccentric¬ 
ity 0.833 and AOP and true anomaly are 0°. The capture 
orbit around Moon is 21,600 x 42,338 km and AOP 0°. 
One of the solutions obtained by present method with GA 
gives AV of the order of 296.5 m/s and time of flight is 
100 days. In order to start from 200 km circular LEO to 
EPO (200 x 65,947 km), AV is 2.755 km/s. So total AV 
is 3.052 km/s. If we start from GTO (200 x 35,840 km) to 
EPO (200 x 65,947 km), AV is 300 m/s and so the total AV 
is 596.5 m/s. 
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Table 2 Comparison between 
Hohmann transfer, invariant 
manifold theory (Topputo et al. 
2004) and present method using 
GA 


Method 

AV (m/s) 


AT (days) 


LEO 

GEO 


Hohmann transfer (Topputo et al. 2004) 

3344 

1177 

6.5 

Invariant manifold (Topputo et al. 2004) 

3081 

914 

49 


3085 

918 

119 


3091 

924 

47 

Present method using GA optimization 

3052 

597 

100 


5 Conclusions 

A Weak Stability Boundary (WSB) trajectory to Moon is 
designed in two stages. First, a highly elliptical geocen¬ 
tric (HEG) orbit is propagated forward in time so that its 
perigee increases to Earth-Moon distance. Second, a highly 
elliptical selenocentric orbit (Lunar capture orbit) is propa¬ 
gated backward in time till it starts moving towards Earth. 
The two trajectories are patched using Fixed Time of Ar¬ 
rival Targeting method to obtain a WSB transfer trajectory. 
It is observed that the Lunar capture trajectories obtained by 
back-propagation is grouped in the phase space with respect 
to time of flight. It is also observed that initial phase angle 
plays an important role in getting captured at Moon. HEG 
orbits are studied in Sun-Earth-Moon-spacecraft R4BP. Tra¬ 
jectories which encounter Moon at 1st perigee passage are 
represented in the phase space. Typical trajectories are also 
shown. Fly-by moon on the way to apogee helps to reduce 
time of flight. Initial lunar phase angle is adjusted to ob¬ 
tain fly-by. Those trajectories which arrive at moon are also 
represented in phase space. Dependence of time of flight 
on angular position of Moon during fly-by is also investi¬ 
gated. Genetic Algorithm is used to find patching points be¬ 
tween HEG and lunar capture orbit. This study will be help¬ 
ful for mission designers. The phase space diagrams with 
color code on time of capture enables selection of departure 
and arrival orbits and an approximation of total flight dura¬ 
tion can be made without actually constructing the complete 
trajectory. 
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Appendix: Equations of motion 

Restricted three-body problem (R3BP) and restricted four- 
body problem (R4BP) are the base models for construc¬ 


tion of weak stability boundary transfer trajectories. In these 
models the orbits of Earth-Moon and Sun-Earth are consid¬ 
ered to be circular and coplanar. But in actual solar system, 
Moon’s orbit eccentricity is 0.0055, Earth’s orbit eccentric¬ 
ity is 0.017 and the Moon’s orbit is inclined to Earth’s or¬ 
bit by 5°. These values are low and so R3BP and R4BP 
can be considered as a good starting model. An orbit which 
becomes a real mission typically is obtained first in such 
an approximate system and then later refined through more 
precise models which include effects such as out-of-plane 
motion, eccentricity, other planet’s gravity, solar radiation 
pressure etc. However, tremendous insight is gained by con¬ 
sidering a simpler model which reveals the essence of the 
transfer dynamics. 

A.l Circular restricted three-body problem (CR3BP) 
(Szebehely 1967) 

It is assumed that two massive bodies called primaries m\ 
(more massive) and m 2 move in circular orbit about their 
center of mass (COM). The third body P is very small to af¬ 
fect the motion of other two bodies. Mass ratio of the system 
is given by 


fju — . y±-rj 

m 1 + m2 

The masses of primaries m\ and m 2 in this non-dimension- 
alized system of units are, respectively, 

/xi = 1 —/x and /X 2 =/x, where /x e [0,1/2]. (15) 

1 unit distance = distance between m\ and m 2 . 

Unit of time is such that the orbital period of m\ and m 2 
about their COM is 2n. 

Let X-Y-Z be an inertial frame with origin at mi-m 2 
center of mass, where X-Y plane is the orbital plane of the 
primaries. Consider the set of axes x,y such that the x-axis 
lies along the line from m\ to m 2 with y-axis perpendicu¬ 
lar to it, completing a right-handed co-ordinate system. The 
x-y frame rotates with respect to the X-Y inertial frame 
with an angular velocity equal to the mean motion, n of 
either mass. We refer to this co-ordinate frame as rotating 
frame. 
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Assume that the two frames coincide at t = 0. Let (X , Y, 
Z) and (x, y, z) be the position of P in the inertial and ro¬ 
tating frames, respectively. In normalized units, we have the 
following transformation of the particle’s position between 
two frames (Koon et al. 2006) 

( X\ /cos t —sin t 0\ /x\ 

Y 1 = I sin t cos t 0 ] I y ) . (16) 

z/ V 0 01/ \z) 


Differentiating gives the transformation of velocity compo¬ 
nents from the rotating to the inertial frame: 


(x\ 

Y 

W 


( COS t 

sin t 
0 


— sin £ 0 \ (x — y\ 

cos t 0 I I y + v I 

0 1A z ) 


(17) 


The Euler-Lagrange equations are given by 

d dL dL _ q 

dt 3 q l 3 q l 

L(x,y,z,x,y,z) = i((i - y) 2 + (j +x) 2 + z 2 ) 

- U(x,y,z), 

y- (x - y) - (y + x) + U x = 0, 
dt 

d . 

— (y + x) + (x - y) + U y = 0, 
dt 

On simplification, we have 


(25) 

(26) 

(27) 

(28) 
(29) 


The larger primary mi, is located at (—/z2,0,0), and the 
smaller primary m 2 , at (/zi,0,0). This is also true in the in¬ 
ertial frame when t = 0. At general time t, the positions of 
mi and m 2 in inertial frame are respectively, 

(Ai,Ti,Zi) = (—/Z 2 cos sin L 0), (18) 

(X 2 , Y 2 , Z 2 ) = (/Ji 1 cos/,/zi sin^O). (19) 

The gravitational potential which the particle P experiences 
due to mi and m 2 (in normalized units) is 

tt ^ 1 /om 

U = --Ml M2- (20) 

r 1 r 2 2 

where 

r\ = (v + /z 2 ) 2 + y 2 +z 2 and r\ = (x - ti\) 2 + y 2 + z 2 . 

( 21 ) 

Consider Euler-Lagrange equations where the mechanical 
system is described by generalized co-ordinates (q 1 ,..., q n ). 


x - 2y = -U x . 

(30) 

y + 2x = —U y , 

(31) 

Ml : 

II 

1 

(32) 

U (x,y,z) = -i(* 2 +y 2 ) + U(x,y,z), 

(33) 


where U is the effective potential and the subscripts de¬ 
note its partial derivatives. Since the equations of motion of 
CR3BP are Hamiltonian and independent of time, they have 
an energy integral of motion. 

E(x,y,z,x,y,z) = ^(x 2 + y 2 + z 2 ) + U(x,y,z). (34) 

Jacobi integral is given by 

7(x,y,z,x,y,z) = -2 E = -(i 2 + y 2 +i 2 ) 

-2 U(x,y,z), (35) 

+ + “2 U(x,y,z)) =0. (36) 


d dL dL _ Q 

dt 3 q l 3 q l 

Lagrangian L = KE — PE 
Lagrangian in Inertial Frame: 


This is the most important integral for determining the mo¬ 
tion of the particle. 

If (p(t) = (x(t),y(t),x(t),y(t)), t g M 1 , represent a gen¬ 
eral solution for (1), then J(cp) = C, where C is a real con¬ 
stant. The set 


C(X, Y, Z, X, Y, Z, t ) = I(X 2 + Y 2 + z 2 ) 

-U(X,Y,Z,t). (23) 

Lagrangian in Rotating Frame: (time independent) 

L(x,y,z,x,y,z) = i((i - y) 2 + (y + x) 2 + z 2 ) 

-U(x,y,z ) (24) 


J(C) = {(xj,ij)eR 4 | J(x,y,x,y) = C], (37) 

For a given value of C, represent a 3D energy surface in the 
4D phase space. 

There are no other integrals constraining the motion of 
the particle, making the CR3BP a non-integrable problem. 

Lagrange points There are five critical points of the effec¬ 
tive potential—three points along the v-axis and two sym- 
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metric points off the v-axis. These points are the x-y lo¬ 
cations of the equilibrium points for a particle in the rotat¬ 
ing frame, i.e., a particle placed here at rest with respect to 
m\ and m 2 (zero initial velocity), will stay at rest for all 
time (zero acceleration). These points are labelled as Li, 
i = 1,... ,5. L\, L 2 and L 3 lie on the line joining the two 
primaries while L 4 and L 5 form an equilateral triangle with 
the primaries. Let Ei be the energy of a particle at rest at Li, 
then £5 = £4 > £3 > E 2 > E\. Thus, L\ is the location of 
the lowest energy equilibrium point and L 4 and L 5 are the 
highest energy equilibrium points. 

A.2 Restricted four-body problem 

A.2.1 Equations of motion in Sun-Earth rotating 
co-ordinate system (Koon et al. 2006) 

In this model we suppose that the Sun and Earth are revolv¬ 
ing in circular orbits around their barycentre and the Moon 
is moving in a circular orbit around the center of Earth. The 
orbits of all four bodies are in the same plane. Let [i, 1 — [i 
and my[ be the masses of Earth, Sun and Moon, respectively. 
Let the distance between the Sun and Earth be unity. The dis¬ 
tance from Earth to Moon is In rotating co-ordinates wrt 
Sun-Earth system the positions of Sun and Earth are fixed at 
and (1 — /i,0), respectively. The angular velocity of 
the Moon in these synodic co-ordinates be com and the phase 
of the Moon at t = 0 be $mo- In rotating frame defined above 
and using non-dimensional units the equations of motion are 

x = x + 2y - cs(x + /ie) - ce(x - /is) ~ cm(x ~ *m), 

(38) 

y = y — 2x — (ce + cm + cs)y + GvOM, (39) 


where c/ = for i = E,M,S, and as = and 

r i a s 


rs = y/(x + Me) 2 + y 2 , 

(40) 

rE = y/(x- Us) 2 + y 2 , 

(41) 

ru—fx- xu) 2 + O - 3'm) 2 - 

(42) 

with 



[is = 1 - [1, [IE —[l, XM = aM COS(^m), 
yM = am sin(0M), 6m = oiMt + 0mo- 

The values of the parameters are [i = 3.0359 x 10 -6 , mM = 
3.733998734625702 x 1(T 8 , a M = 2.573565073532068 x 
1(T 3 and com = 12.36886949284508. 


A.3 Fixed time of arrival targeting (Yamakawa 1992) 

Suppose at initial epoch (to), position and velocity (r(to), 
v(to)) are given. 

Transition matrix from to to tf is 

*^={Z %)• (43) 

where @ij are 3 x 3 partitioned transition matrix. 

(££!)• (44) 

where dV 0 is velocity correction at to and dr and dv are 
position and velocity deviations, respectively. 

The objective is to cancel dr(tf) + , 

< &n(tf,to)dr(to) + <Pn(tf,to){dv(to) -\-dVo) =0 

which yields 

&i 2 (tf,to)dVo = —&n(tf,to)dr(to) 

~ <P\ 2 (tf,to)dv(to), 

dV 0 = —<Pi 2 (tf,to)~ l [<Pn(tf,to)<Pi 2 (tf,to)]x(to). 

When dVo is not performed at to, positional deviation at tf 
is 

dr(tf)~ = [@n(tf,to)d>i 2 (tf,to)]x(to). 

Thus, velocity correction at to is given as 

dV 0 = -0n(tf,to)~ l dr(tf)~. (46) 

The time required to travel from starting point to end point 
is optimized using Nelder-Mead simplex method Lagarias 
et al. (1998) to minimize the total impulse required. 
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